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We report internally contracted relativistic multireference configuration interaction (ic-MRCI), complete ac¬ 
tive space second-order perturbation (CASPT2), and strongly contracted /i-electron valence state perturbation 
theory (NEVPT2) on the basis of the four-component Dirac Hamiltonian, enabling accurate simulations of rela¬ 
tivistic, quasi-degenerate electronic structure of molecules containing transition-metal and heavy elements. Our 
derivation and implementation of ic-MRCI and CASPT2 are based on an automatic code generator that trans¬ 
lates second-quantized ansatze to tensor-based equations, and to efficient computer code. NEVPT2 is derived 
and implemented manually. The rovibrational transition energies and absorption spectra of HI and T1H are 
presented to demonstrate the accuracy of these methods. 


I. INTRODUCTION 


There are continued interests in accurate modeling of gas- 
phase thermochemistry and dynamics that involve transition- 
metal and heavier elements, where relativistic effects play 
an important role. For instance, scientists at the US 
Air Force recently performed an experiment, aiming to 
use chemi-ionization involving lanthanide atoms to alter 
the electron density in the ionosphere for radio-frequency 
communication, 12 for which accurate simulations could help 
analyze the experimental observation. Another example is the 
reaction of FeO + with a hydrogen molecule, a model reac¬ 
tion system for the so-called two-state reactivity,® of which 
accurate modeling still remains a challenge!- 4 The spin bar¬ 
riers in the two-state reactivity mechanism are also ubiqui¬ 
tous in organometallic chemistry.® Understanding these prob¬ 
lems requires accurate description of strongly relativistic, 
quasi-degenerate electronic structure. There have been, how¬ 
ever, only a handful of theory developments to address this 
challenge.®® 

As a first step toward realizing predictive simulations of 
such processes, we develop in this work novel computa¬ 
tional tools that combine the four-component relativistic Dirac 
formalistrP- 1 and internally contracted multireference electron 
correlation methods. Our approach is based on the four- 
component Dirac equation for electrons. 
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where a and f3 are Dirac’s matrices, g(i, j) is a two-electron 
operator, and c is the speed of light. Z A is a charge of a 
nucleus A (note, however, that we use finite-nucleus models 
in practice). Hereafter atomic units are used unless other¬ 
wise stated. In this work, we use the full Breit operator for 
electron-electron interactions, i.e.. 
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The reader may consult Ref. Jl2]for details on integral evalu¬ 
ation associated with this operator over Gaussian basis func¬ 
tions. We first perform complete active space self-consistent 


field (CASSCF) calculations using this Hamiltonian, '31131 j n 
which orbitals are optimized using the minimax principle, and 
project out the space spanned by the ‘negative-energy’ or¬ 
bitals, a procedure called no-pair projection. 11 An efficient 
Dirac-CASSCF algorithm that we have developed can be 
found in Refs. [14] and [13] After the no-pair projection pro¬ 
cedure, the Hamiltonian in the second quantization becomes 


#NP - ^syExy + ^ 

xy xyzw 
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where x, y, z, and w label any electronic molecular spin 
orbitals (MO), and h xy and v xy , tV , v are the (complex-valued) 
Hamiltonian matrix elements in the MO basis in chemists’ no¬ 
tation. E xy and E xyjw are operators defined as 

E xy ti^ciy, (4a) 

Exy,zw — a\a\a w a r (4b) 

Since the MO Hamiltonian [Eq. Q] is isomorphic to the 
non-relativistic counterpart (and all the eigenstates are min¬ 
ima in the parameter space after the no-pair projection pro¬ 
cedure), standard electron-correlation methods, such as inter¬ 
nally contracted multireference configuration interaction (ic- 
MRCI),®® 9 can be used in conjunction with this Hamilto¬ 
nian. We note in passing that, even though our numerical re¬ 
sults are based on the four-component formalism [Eq. 0]. 
the multireference theory and programs developed in this 
work are equally applicable to any two-component relativistic 
Hamiltonians.®® 

In the non-relativistic framework, the ic-MRCI method has 
been pioneered by Werner and co-workers.®® The ability 
of ic-MRCI to accurately and consistently describe the poten¬ 
tial energy surfaces of small-molecule reactions has been the 
key to understanding many of the gas-phase reactions studied 
in the past decades (for instance, see Refs. [2314251 . Very re¬ 
cently ic-MRCI has been extended to incorporate density ma¬ 
trix renormalization group reference functions with more than 
20 orbitals in the active space by Saitow et al. 19 There are also 
parallel implementations of uncontracted MRCI,® though its 
computational cost is generally higher than that of ic-MRCI. 

Another class of popular multireference approaches in non- 
relativistic theory is based on perturbation theory. Among 
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others the complete active space second-order perturbation 
(CASPT2) methocP® is an internally contracted, multiref¬ 
erence generalization of the standard Mpller-Plesset pertur¬ 
bation theory and has been applied to a wide variety of chem¬ 
ical problems. 1 ’ 0 The ^-electron valence state perturbation the¬ 
ory (NEVPT2) 3132 proposed by Angeli et al. (especially 
its strongly correlated variant) uses a different zeroth-order 
Hamiltonian and has desirable properties such as strict size ex- 
tensivity and numerical robustness against so-called intruder- 
state problems. 

Here we report the theory and algorithms for relativistic ic- 
MRCI, CASPT2, and NEVPT2 based on the four-component 
Dirac Hamiltonians. This work realizes relativistic ic-MRCI 
and NEVPT2 for the first time, whereas CASPT2 has been 
reported in the past by Abe et al. 7 and by Kim et al.® The 
implementations of ic-MRCI and CASPT2 are facilitated by 
an automatic code generator, smith3. 33 ® The smith3 program 
was previously used to derive and implement nuclear en¬ 
ergy gradients for fully internally contracted CASPT2 33 and 
has been extended in this work to incorporate equations with 
spin orbitals in complex arithmetic. Note that the automatic 
code generation approach has been used for relativistic single¬ 
reference coupled-cluster methods by Hirata et al.® and by 
Nataraj et al. 36 The generated code and the code generator are 
both publicly available EC 37 The NEVPT2 code is manually 
implemented. In the following we sketch the outline of the 
theories and implementations. 

II. THEORY 

A. Relativistic MRCI with internal contraction 

Our ic-MRCI implementation uses fully internally con¬ 
tracted basis functions, which are similar to those used in the 
CASPT2 theory by Roos and co-workers.® The correlated 
wave functions are parameterized as 

|T> = 7ref|<D> ref ) + 2 T n E n \® ref), (5) 

n 

in which 7”s are the unknown amplitudes to be determined, Q 
denotes excitation manifolds in ic-MRCI, and E<_> are associ¬ 
ated excitation operators: 

Co — | E a ibj, E„r,hi, E ar jy S , E mr] , 

Eri'Sji E ar s t , E r i St , E a [ rs I. (6) 

Hereafter i and j label closed orbitals, r, s, and t label active 
orbitals, and a and b label virtual orbitals. Note that, because 
spin orbitals are used, E m ,rs and E ru ,n that are distinguished in 
non-relativistic theories generate identical sets of excited con¬ 
figurations. The Kramers symmetry is not utilized in our ic- 
MRCI implementation except for integral compression. |O re f) 
is a relativistic multi-determinant reference function, 

|G> ref >= J] C" + ’"-|/" + '"->, (7) 

n+ +ri-=n 


where n+ and n are the numbers of electrons that belong to 
Kramers + and - spin orbitals, and n is the total number of 
active electrons! 13 ! 14 ! 

In the ic-MRCI method, the Dirac Hamiltonian is diagonal¬ 
ized in the space spanned by the parameters in Eq. ®,i.e., 

E = min [<'T|H np |'P>] , (8) 

under a normalization constraint. The following cr and n vec¬ 
tors are computed from each trial vector <///> in the same basis. 


(o' p)n = (O re f |£^H NP t/q>). 
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(cp)ref = ( < Sref l#NP 
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(xp)n = (‘I'ref Eq \^p). 
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Note that we eliminate five-particle reduced density matrices 
from the equations by means of a well-known commutator 
trick, i.e., (using T<> = TqEo) 

( Oi-ef | H NP fn | <Dref ) 

= (4 ) ref|£ , Q<[HNP, Tnll^ref) + (‘IWIE’q, TnlOref^ref, (10) 

where Q and Q' belong to the same excitation class in Eq. ([6}. 
A Hamiltonian matrix is then constructed within the subspace 
spanned by the trial vectors,® 

HpQ—TpCTQ, SpQ=Tp7tQ, (11) 

and diagonalized to obtain the coefficients (cp) that constitute 
an optimal linear combination of the trial vectors: 

^ h pq c Q = E ^ SpqCq. (12) 

Q Q 

Using these quantities, the residual vectors are 

R = ^ Cp[(Tp - E7T p ] , (13) 

p 

from which we generate a new set of trial vectors (see below). 

The working equations [Eqs. (|9a|)-([T0|)] for cr-vector for¬ 
mation can be expressed in terms of reduced density matrices; 
therefore, it is essentially identical to the non-relativistic coun¬ 
terpart except for spin symmetry in the latter. The explicit 
formulas consist of ca. 750 tasks, most of which are tensor 
contractions. They can be found in supporting information. 39 
The equations were implemented into efficient computer code 
using the automatic code generator smith3. 33 ® First, smith3 
performs Wick’s theorem to convert second-quantized expres¬ 
sions to a list of diagrams represented by tensors and their 
contractions. Next it factorizes the diagrams to a tree of binary 
tensor contractions. Finally the tree is translated to computer 
code that is compiled and linked to the bagel package. 37 See 
Refs. I4QM421 for further information on automatic code gener¬ 
ation. 
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At the end of each ic-MRCI calculation, the Davidson cor¬ 
rection is added to the total energy to approximately account 
for size-extensivity errors.®] The correction is 


AE+q - 


T 2 \ 
1 ref 


T 2 

1 ref 


(14) 


where T ccf is the weight of the reference configuration in the 
correlated wave function [see Eq. <|5j], and E con is the corre¬ 
lation energy from ic-MRCI calculations. 


B. Relativistic CASPT2 and NEVPT2 


The second-order perturbation methods, CASPT2 and 
NEVPT2, are defined as minimization of the so-called Hyller- 
aas functional, 

E = min [<'P (1) |// (0) - £ ,(0) |'F C1) > + 2^<'P a) |H NP |<D ref )]. (15) 

In CASPT2, the zeroth-order Hamiltonian H (0) is chosen to be 
a projected Fock operator 

H (0) = PfP + QfQ, (16) 


where P is a projector to the reference configuration and Q 
is its orthogonal compliment. The first-order wave function 
I T (| J is parameterized as in Eq. 0 The minimization is per¬ 
formed by solving a set of linear equations using a subspace 
algorithm. The construction of residual vectors, 

Rn = 2 [<Q|H (0) - £ ( 0 Vp> + <Q|H N p|<l>ref)], (17) 


is akin to (but simpler than) that in ic-MRCI. Here we used 
(Q| = (d> rc ijT(,. For details on the relativistic CASPT2 equa¬ 
tions, see earlier reports by Abe et al. 7 and Kim et a 1 P®] 

In NEVPT2, the zeroth-order Hamiltonian is defined using 
Dyall’s Hamiltoniarf^as 

H (0) = PH np P + 2 |® U >£ W <®J, (18) 


where a> is the excitation class in Eq. 0 and ct> w is defined as 


l®«> = 


PtoH NP |<E) re f) _ 

V( < l ) ref I^NpP w^NP^ref) 


(19) 


P ( „ is a projector onto a>, and the denominator accounts for 
normalization. E l:> that appears in Eq. © is 

Eo, = <® ( JHnpI® w >. (20) 


The wave function is parameterized using the so-called strong 
contraction scheme, i.e., 

IT) = T rcf |O rcf ) + ^ 7J®„>. (21) 

at 

Since // ,0i of NEVPT2 does not include off-diagonal cou¬ 
plings between different ai, the equations can be solved with¬ 
out iterative procedures. The working equations for relativis¬ 
tic NEVPT2 can be obtained by dropping the factors of 2 that 
stem from spin summations in the non-relativistic equations 
in Ref. [32] The explicit formulas are provided in supporting 
information. 39 


C. Wave function updates in ic-MRCI and CASPT2 

Internally contracted basis functions (£n|O re f)) are not or¬ 
thogonal with each other and sometimes linearly dependent;®] 
therefore, one has to take into account the overlap matrix 
when updating the amplitudes. The generation of trial vectors 
is performed as the following. Let us consider as an example 
the amplitudes associated with E ar \, s . In this case, the overlap 
and (approximate) diagonal Hamiltonian matrix elements, S 
and F, respectively, are 

Srs,r's' ~ (TVcdEVr'.vv'|4^ rc | ), (22a) 

Ers.r's 1 ~ ^ ' .(^rci'\EJt' [4\ c l).//r' 5 (22b) 

11' 


where E rr , ss p tV = a r E ss , tt :a r i. We calculate while pro¬ 
jecting out the linearly dependent part so that (S -1 ^ 2 ) 1 SS 1/2 is 
a unit matrix (the eigenvalues that are smaller than 1.0 x 10 8 
are discarded), which is then used to form 

F = (S -1/2 ) t FS _1/2 . (23) 

Next F is diagonalized to yield a transformation matrix U, 

F = U/lU t , (24) 

with a diagonal matrix A. Defining X = U^S -1 ^ 2 , we arrive 
at the formula for generating new trial vectors from residual 
vectors: 


(>j / p +1 )ar,bs - ^ y. 


Par’ Jis' A1),- 


E {0) - A d - e a - e h 


X 




(25) 


where e a is an orbital energy (i.e., e a = f aa ) and D labels the 
eigenvalues in Eq. ( |24| ), the number of which is equal to or 
smaller than the numbers of rows and columns of the overlap 
matrix [Eq. ©]• This formula implies that in ic-MRCI up¬ 
dates the inverse of Hnp - E is approximated by that of the 
diagonal part of the CASPT2 equation. 27 


D. Computation of rovibrational spectra 


Rovibrational energy levels of diatomic molecules in their 
X states can be calculated by solving an effective one¬ 
dimensional Schrodinger equation (in this section we avoid 
use of atomic units for clarity), 


n 2 d 2 

2/r dr 2 


+ V(r ) + -— ^J(J + 1) 
2/jr- 


TV/ (r) = Ey j'Vy j (r ), 

(26) 


in which v and J are the vibrational and rotational quantum 
numbers, respectively, and p is the reduced mass. The third 
term of the Hamiltonian accounts for the Coriolis coupling. 
The rotation-vibration coupling is, therefore, variationally in¬ 
cluded in the calculations. 
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TABLE I. Root-mean-square deviations of the rovibrational transi¬ 
tion energies of H 127 I and 205 T1H in cm' 1 computed by the four- 
component methods. The HITRAN database^ and experimental 
dad® were used as references. 




CASSCF 

CASPT2 

NEVPT2 

MRCI+Q 

Origin 

HI 

v — 0 —> 

1 

120 

36 

21 

8 

2230 

T 

O 

II 

2 

245 

73 

43 

15 

4379 

T 

O 

II 

3 

378 

117 

68 

24 

6448 

T 

O 

II 

4 

519 

163 

96 

34 

8435 

T1H 

v = 0 —> 

1 

92 

34 

47 

17 

1345 

v = 1 —> 

2 

93 

33 

46 

15 

1300 

v = 2 —7 

3 

92 

33 

46 

13 

1255 



Bond length / A 


FIG. 1. Potential energy curves of HI computed by four-component 
CASSCF, CASPT2, NEVPT2, and ic-MRCI+Q.^ The experimental 
bond length and dissociation energy are 1.609 A and 3.20 eV. re¬ 
spectively. 


The line intensity Iy associated with the transition energy v 
can be computed ad^S 


/* = 


(2 J f + 1) 

8 ttcQv 2 


My,. 




-Ei/kT 


(1 


-hcv/kT\ 


h 


(27) 


in which E ,• is the energy of the initial state and k is the Boltz¬ 
mann constant. The partition function Q at a temperature T is 
evaluated using 


Q = Yj(2J, + l)e- E ‘ /kT , (28) 

/ 

where / runs over rovibrational states. We used T = 296 K. 
The quantum numbers of initial (final) states are labeled by 
v, and Jj (vj and Jf). Using the rovibrational wave functions 
( V P v .j. and j ) and the dipole-moment function Mir), the 

Einstein coefficient J[ Vj j.^ Vf Jf is 

87r 2 v 3 S^ j, . ,2 

- 3w 3 fl2J|+ 'i l<^IM W |T^)| . (29) 

where eo is the vacuum permittivity and S j t j. is the Honl- 
London factor,® which is max(7,, Jf) for the electronic 
ground states of HI and T1H. 


III. NUMERICAL RESULTS 

First, to benchmark the accuracy, we applied four- 
component CASSCF, CASPT2, NEVPT2, and ic-MRCI+Q to 
an HI molecule, for which there are reliable experimental ref¬ 
erence data. 47 Uncontracted Dyall’s cv3Z® and uncontracted 
cc-pVT2T° basis sets were used for I and H, respectively. 
Gaussian-type nuclear charge distributions were used. 51 The 
4,y, 4p, 4 d, 5s, and 5p electrons of I and the 1 ,v electron of 
H were correlated (i.e., 26 correlated electrons; 28 electrons 
were frozen), among which 5,v, 5p of I and I s of H were 
treated in the active space. In correlated calculations, virtual 
orbitals were truncated at 55 E^. The total number of corre¬ 
lated spin orbitals was 206. The computed potential energy 


curves relative to their minima are shown in Fig. |T] The equi¬ 
librium bond lengths obtained by CASPT2, NEVPT2, and 
ic-MRCI+Q were 1.608, 1.609, and 1.606 A, respectively, 
which are in good agreement with the experimental value 
(1.609 A)P 2 The dissociation energies D e were estimated via 
extrapolation to be 3.0, 3.0, and 3.1 eV, respectively. The ex¬ 
perimental value is 3.20 eV. 52 

We then simulated the absorption spectra based on these 
potential energy curves interpolated by five-point piece-wise 
polynomials. Dipole moments were computed at each point 
as electric-field derivatives [ M(r ) = dE(r)/d& : where S z is 
an external electric field along the molecular axis] using finite 
difference formulas. The Level 8.2 progranl® was used to 
solve the radial Schrodinger equation [Eq. ( |26[ i | and to evalu¬ 
ate [Eq. (|29]i |. The partition function and absorp¬ 

tion spectra were computed using a program of Yorke et al.® 
The computed spectra for the fundamental, overtone, and sec¬ 
ond overtone transitions are presented in Fig. [2] in which the 
HITRAN reference spectra 47 are also shown. Overall, the line 
positions were accurately reproduced by ic-MRCI+Q within 
0.5 % (8 cm -1 for the fundamental transitions and 34 cm -1 for 
the third overtone transitions), attesting to the consistent ac¬ 
curacy of ic-MRCI+Q throughout potential energy surfaces; 
The line intensity of the overtone and second overtones agreed 
well. Our results overestimated the intensity of the fundamen¬ 
tal transitions, which is mainly because the intensity is largely 
suppressed by the almost flat dipole-moment curve around the 
equilibrium geometry; therefore, it is highly sensitive to the 
accuracy of the computed dipole moments.® The errors in 
the line positions computed by CASPT2 and NEVPT2 were 
found three or four times larger than those by ic-MRCI+Q. 

Next, we calculated the potential energy curve of T1H us¬ 
ing CASSCF, CASPT2, NEVPT2, and ic-MRCI+Q. The elec¬ 
tronic structure of T1H around the equilibrium geometry has 
been studied by many authors .Ss M H We used uncontracted 
Dyall’s cv 3 /® and uncontracted cc-p VT2P basis sets for T1 
and H, respectively, in conjunction with Gaussian-type nu¬ 
clear charge distributions.® The full-valence active space (4 
electrons in the 6s and 6p orbitals of T1 and the Is orbital 
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FIG. 2. Simulated rovibrational absorption spectra of H I27 I at 296K using four-component CASSCF, CASPT2, and ic-MRCI+Q. The bottom 
panels are the observed lines from the HITRAN database (hyperfine-split lines are averaged for comparison). 



Bond length / A 

FIG. 3. Potential energy curves of T1H computed by four-component 
CASSCF. CASPT2, NEVPT2, and ic-MRCI+Q. The experimental 
bond length and dissociation energy are 1.872 A and 2.06 eV, re¬ 
spectively. 


of H) was used. The 5s, 5 p, 4 /, 5 d, 6s, and 6 p electrons of 
T1 and the Is electron of H were correlated (i.e., 36 corre¬ 
lated electrons). The virtual orbitals were again truncated at 
55 £h. resulting in 248 correlated spin orbitals. The potential 
energy curves of T1H computed by four-component CASSCF, 
CASPT2, NEVPT2, and ic-MRCI+Q are shown in Fig. [3] The 
dissociation energy D e from ic-MRCI+Q (2.00 eV) was in 
excellent agreement with the experimental value (2.06 eV)J^ 


while CASPT2 underestimated it by 0.2 eV (1.84 eV). The 
equilibrium bond length (1.872 A) was also accurately re¬ 
produced by ic-MRCI+Q (1.872 A). Those by CASPT2 and 
NEVPT2 were 1.870 and 1.885 A, respectively. NEVPT2 was 
found less accurate than CASPT2 for this molecule, and its 
accuracy deteriorated as the bond is stretched. 

The absorption spectra of T1H were likewise computed us¬ 
ing the energies at 20 grid points between 1.3 A and 6.0 A. 
The computed spectra are presented in Fig. [4] The exper¬ 
imental line intensity was not found in the literature. The 
mean-root-square errors in the computed rovibrational tran¬ 
sition energies are also listed in Table. [T} in which the ex¬ 
perimental results from Ref. [48] are used as reference val¬ 
ues. The errors in the transition energies were around 35, 45, 
and 15 cm 1 forCASPT2, NEVPT2, and ic-MRCI+Q. Apart 
from the shift, the line positions computed by ic-MRCI+Q 
agree perfectly with the experimental results. The remain¬ 
ing errors include incomplete treatment of dynamical corre¬ 
lation in the ic-MRCI+Q model, the effects of the higher- 
order quantum-electrodynamics interactions, and the non- 
Born-Oppenheimer contributions. 

The wall times for one iteration of relativistic CASPT2 and 
ic-MRCI on T1H were roughly 2 and 80 minutes using two 
Xeon E5-2650 CPUs (2.0 GHz, 8 cores each) on a single 
node. The wall time for non-relativistic ic-MRCI per iteration 
is about 16 seconds; therefore, relativistic ic-MRCI is roughly 
300 times more expensive than the non-relativistic counter- 
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FIG. 4. Simulated rovibrational absorption spectra of 205 TIH at 296K 
using four-component CASSCF, CASPT2, and ic-MRCI+Q. Dotted 
lines in the bottom panel are the experimental line positions taken 
from Ref. |48] superimposed by shifted ic-MRCI+Q spectra. 


part. A factor of 2 6 = 64 stems from the fact that relativistic 
ic-MRCI does not use spin symmetry. An additional factor 
of 3 should be ascribed to matrix multiplication in complex 
arithmetic that is three times as expensive as that in real arith¬ 
metic. The rest is due to other factors such as caching and 
optimized libraries. 


IV. CONCLUSIONS 

In summary, we have developed four-component relativis¬ 
tic ic-MRCI, CASPT2, and NEVPT2 based on the Dirac 
Hamiltonian and full internal contraction. The relativistic ic- 


MRCI and CASPT2 programs have been implemented using 
automatic code generation. The programs are interfaced to 
the open-source bagel package. 37 The code generator smith3 
is also publicly available.® The accuracy of these methods 
has been presented by computing the entire potential energy 
curves of HI and T1H and directly comparing calculated rovi¬ 
brational transition energies with the experimental data. It 
has been shown that ic-MRCI+Q can reproduce experimen¬ 
tal transition energies with 0.5 % and 1 % accuracy for HI 
and T1H, respectively, up to high-lying rovibrational transi¬ 
tions using uncontracted triple-^ basis sets without any cor¬ 
rections or extrapolations. 

Currently the size of ic-MRCI and CASPT2 calculations is 
limited by the memory requirement for two-electron MO in¬ 
tegrals that are stored in core, which is somewhat problematic 
especially because uncontracted one-electron basis functions 
(with energy cut-offs) have to be used for heavy elements. 
Furthermore, wall times for multi-state ic-MRCI calculations 
scale cubicly with respect to the number of states included in 
the calculation, which become prohibitively long when sev¬ 
eral states are included in the calculations. To address these 
problems, the parallelization of the programs based on the 
tiledarray library of Calvin and Valeev® is under develop¬ 
ment in our group. Our relativistic NEVPT2 code does not 
store 4-index intermediates and is heavily parallelized (to be 
presented elsewhere); therefore, it is ready for use in chemical 
applications. 

SUPPORTING INFORMATION 

The working equations for relativistic NEVPT2 and the 
rovibrational transition energies and absorption spectra of 
HI and T1H can be found in supporting information. The 
computer-generated ic-MRCI equations are also included. 
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